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—.5.  Contours  are  cumulative  probability  that  the  intercept  slope  pair  is  contained 
within  the  region. 
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1.  Introduction 


There  are  two  important  problems  in  monitoring  nuclear  tests  using  seismic  and  other 
means,  either  in  a  proliferation  or  treaty  monitoring  environment.  The  first  is  calibrating 
an  assumed  relation  between  the  yield  of  a  nuclear  test  and  a  vector  of  observed  magnitudes. 
For  example,  we  may  measure  body  wave  magnitudes,  mb  and  Lg  for  a  number  of  nuclear 
tests  and  seek  to  calibrate  some  function  to  the  observed  data  which  will  allow  prediction 
of  yield  from  the  magnitude  vector.  The  prediction  of  yield  and  its  associated  uncertainty 
from  observations  on  a  future  magnitude  vector  is  the  second  important  problem  of  interest 
considered  in  this  report. 

To  illustrate  the  nature  of  the  first  problem,  consider  a  set  of  16  explosions  from 
the  Russian  site  at  Semipalatinsk  for  which  yields  have  been  published  in  Bocharov  et  al 
(1989)  and  Vergino  ( 1989).  The  data,  listed  in  Table  1  of  Section  5,  are  plotted  in  Figure  1 
where  it  is  fairly  clear  that  a  linear  model  relating  both  magnitudes  to  the  logarithm  of  the 
announced  yield  would  not  be  unreasonable.  The  only  apparent  deviation  from  linearity 
is  the  second  event  where  the  mb  value  seems  out  of  line.  Fitting  linear  regressions  of  the 
two  magnitudes  on  log-yield  seems  to  be  a  natural  way  of  estimating  a  calibration  relation. 

Note  that  all  of  these  events  are  from  the  same  region  and  it  is  unlikely  that  monitoring 
in  a  proliferation  environment  could  yield  such  a  large  calibration  set.  Furthermore,  the 
slopes  and  intercepts  may  vary  considerably  depending  on  the  location  of  the  event  and 
the  location  of  the  network  monitoring  the  event  (see  Heasler  et  al,  1990).  One  would  more 
likely  have  available  a  much  smaller  group  of  say  5  or  6  vector  magnitudes  with  announced 
yields  as  well  as  some  prior  information  about  the  nature  of  the  slopes  and  intercepts. 
The  prior  values  could  based  on  geologic  factors  and  (or)  expert  opinion.  Given  this  more 
limited  situation,  one  could  still  calculate  a  classical  linear  regression  of  the  magnitude 
vector  on  yield  that  would  make  no  use  of  prior  knowledge.  One  could  also  ignore  the 
calibration  data  completely  and  just  use  some  prior  probability  distribution  for  the  slope 
and  intercept  values  based  on  an  uncertainty  derived  from  expert  opinion.  Alternatively, 
there  is  a  third  Bayesian  calibration  approach  that  makes  optimal  use  of  both  the  data 
and  prior  information.  We  will  investigate  all  three  approaches  in  this  report. 

In  the  past,  there  have  been  close-in  hydrodynamic  measurements  called  CORRTEX 
made  on  yield;  these  measurements  are  generally  thought  to  be  more  accurate  than  seismic 
magnitudes.  One  can  either  regard  these  measurements  as  highly  accurate  magnitude 
surrogates  with  unit  slope  and  zero  intercept  or  as  errors  in  variables  measurements  on 
yield.  We  discuss  a  multivariate  approach  using  CORRTEX  as  an  additional  measurement 
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on  yield  in  Section  4.3. 

Yield  estimation  on  the  basis  of  the  calibration  relation  derived  above  is  the  second 
important  problem  we  consider  in  this  report.  We  assume  that  new  yields  corresponding  to 
vector  magnitudes  are  fixed  and  unknown.  Then,  the  predictive  distribution  (see  Atchison 
and  Dunsmore,  1975)  of  the  magnitudes  in  both  the  classical  and  Bayesian  frameworks 
will  be  functions  of  the  univariate  yield.  Inverting  the  predictive  magnitude  distribution 
then  gives  a  confidence  interval  on  the  unknown  future  yield.  An  approach  using  the 
classical  (non- Bayesian)  method  of  Fieller  (1954)  in  the  univariate  case  was  used  in  Rivers 
et  al  (1986)  and  in  Picard  and  Bryson  (1992).  The  extension  of  the  classical  method 
to  the  multivariate  case  is  due  to  Brown  (1982).  For  the  Bayesian  case,  we  choose  to 
regard  the  yields  as  being  fixed  and  unknown  rather  than  as  random  with  uniform  prior 
distributions  as  in  Brown  (1982)  and  Picard  and  Bryson  (1992).  The  use  of  the  predictive 
distribution,  without  integrating  over  yields,  gives  a  classical  confidence  interval  rather 
than  the  Bayesian  interval  which  would  be  characterized  by  the  conditional  distribution  of 
yield  given  the  magnitude  vector  in  the  usual  Bayesian  approach. 

To  summarize  the  approach  taken  here,  we  have  integrated  material  from  Shumway 
and  Der  (1990)  and  Shumway  (1990)  with  the  emphasis  changed  from  monitoring  a  possible 
Threshold  Test  Ban  Treaty  (TTBT)  to  monitoring  the  possible  worldwide  proliferation 
of  nuclear  weapons.  A  new  section  is  included  that  covers  estimating  the  yields  when 
CORRTEX  measurements  are  available  using  maximum  likelihood  and  a  Bayesian  set 
of  priors  for  the  intercept,  slope  vector.  This  allows  yield  estimation  using  only  prior 
information  and  no  data.  To  accomplish  this,  we  first  introduce  the  simple  linear  model 
relating  elements  of  the  magnitude  vector  to  log  yield.  Included  in  this  specification  are 
the  slope  and  intercept  uncertainty,  formulated  in  one  case  using  a  bivariate  normal  prior, 
and  the  magnitude  uncertainty,  formulated  using  the  multivariate  inverted  Wishart  (chi- 
square)  distribution.  We  then  show  how  the  slope  and  intercept  can  be  calibrated  using 
classical  (data  only),  Bayesian  (prior  only  or  prior  combined  with  data)  and  CORRTEX 
observations. 

A  second  Bayesian  approach  is  introduced  that  allows  more  generality  in  formulating 
the  prior  information  on  slopes  and  intercepts  but  assumes  that  the  covariance  matrix 
is  fixed  and  known.  For  yield  estimation,  there  is  a  predictive  approach  using  both  the 
classical  and  Bayesian  assumptions  and  a  Bayesian  approach  using  the  empirical  Bayes 
estimators  for  the  slopes  and  intercepts. 
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2.  Models  for  Magnitude  Data  and  Prior  Information 

Figure  1  implies  that  a  linear  model  relating  log-yield  to  P-vvave  r»*  and  Lg  magnitudes 
is  reasonable,  and  evidence  from  other  research  (see  Ringdal  and  Marshall,  1989)  indicates 
that  other  kinds  of  magnitude  measures  behave  in  a  similar  manner.  Then,  the  linear  model 
relating  the  vector  of  p  magnitudes  from  the  ith  event,  m,  =  (mu , . . . ,  m,p )',  i  —  1 . . . ,  N 
to  the  yields  Wi,  i  =  1, . . . ,  Ar  can  be  written  in  the  form 

m,  =  a  +  bu),  -(-  e,,  (2.1) 

where  a  is  the  pxl  vector  of  unknown  intercepts  and  b  is  the  pxl  vector  of  slopes.  The 
error  vectors  e,  are  assumed  to  be  independent  zero-mean  multivariate  normal  random 
variables  with  common  covariance  matrix  E. 

It  is  convenient  to  collect  the  sample  of  observed  vector  magnitudes  in  the  overall  pxN 
matrix  Y  =  (mi,...,mjv)  and  to  collect  the  error  vectors  in  a  similar  matrix,  denoted 
here  by  V  =  (e1( . . .  ,e/v).  Then,  defining  the  matrices 

X  =  (  1  1  ■"  M  (2.2) 

V  wi  w2  wN  J 

and 

B  —  (a,  b)  (2.3) 

leads  to  writing  (2.1)  as  the  overall  multivariate  linear  model 

Y  =  BX  T  V,  (2.4) 

as  is  given  for  example,  in  Anderson  (1986,  Chapter  8). 

If  there  are  a  limited  number  of  data  points,  it  may  be  important  to  bring  in  prior 
information  about  the  intercept  and  slope  components  of  B  and  the  error  covariance  matrix 
E  .  Suppose,  for  example,  that  the  uncertainty  in  the  covariance  matrix  E  can  be  expressed 
in  terms  of  the  an  inverted  Wishart  distribution  (see  Anderson,  1984,  Section  7.7.1)  with 
parameters  m  and  =  mSo  where  So  is  the  prior  covariance  matrix  corresponding  to  the 
assumed  prior  value  of  E  and  m  is  the  parameter  expressing  the  initial  uncertainty.  The 
value  of  m  corresponds  roughly  to  the  degrees  of  freedom  or  equivalent  sample  size  of  the 
prior  information.  That  is,  how  large  a  sample  would  have  been  required  to  produce, for 
example,  the  uncertainty  of  a  panel-furnished  distribution?  The  prior  distribution  of  the 
covariance  matrix  used  is  of  the  form 
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(2.5) 


7 r(E)  oc  | sr^m+P+2)exp{-^tr  E_1E0} 
where  oc  denotes  proportional  to  and  tr  denotes  the  trace  of  a  matrix. 

In  order  to  illustrate  some  of  the  considerations  involved  in  specifying  the  joint  prior 
distributions  described  above,  consider  the  case  where  we  record  only  one  magnitude  (p=l ) 
and  suppose  that  the  distribution  of  the  initial  variance  is  inverted  Wishart  (chi-square) 
with  m=10  or  m=40  and  <7$  =  .052.  Figure  2  shows  the  probability  densities  of  the 
standard  deviation  a  corresponding  to  the  cases  m=10  and  m=40.  It  is  seen  that  the 
choice  m=10  corresponds  to  a  wider  uncertainty  interval  (roughly  .03  to  .12)  whereas  the 
choice  m=40  narrows  the  interval  (to  roughly  .04-. 07). 

One  must  also  specify  the  uncertainty  of  the  intercept  and  slope  components.  The 
slope  and  intercept  components  in  B  =  (a,  b)  are  assumed  to  be  (conditional  on  E)  nor¬ 
mally  distributed.  The  stacked  vector  consisting  of  the  rows  of  B  is  assumed  to  be  normally 
distributed  with  mean  components  (ao,bo)  and  covariance  matrix  E  0  Dq1  where  0  de¬ 
notes  the  Kronecker  or  direct  product.  The  2x2  matrix  Do  is  chosen  to  reflect  the  joint 
uncertainty  structure  of  a  and  b.  It  corresponds  with  the  joint  uncertainty  structure  of  the 
maximum  likelihood  estimator,  B:  given  above  in  Equation  (3.1)  for  which  the  comparable 
stacked  vector  has  covariance  matrix  E  @C-1 .  For  completeness,  we  give  the  prior  density 
of  the  slope-intercept  matrix  as 

jri(J?|E)  oc  |E| -^exp{-l-tr  E ~\B  -  B0)D0(B  -  B0)'}  (2.6) 

In  Figure  3,  we  see  the  case  corresponding  to  the  intercept  a  and  slope  b  having 
expectations  ao  =  4.4  and  bo  =  9  and  standard  deviations  .2  and  .1.  The  correlation 
was  taken  as  -.5  so  that  the  increase  in  slope  must  necessarily  lead  to  a  decrease  in  the 
intercept.  In  this  case,  the  initial  slope  uncertainty  ranges  roughly  between  .7  and  1.1  with 
the  intercept  ranging  from  3.9  to  4.9.  One  would  expect  that  expert  input  from  panels 
and  studies  would  be  used  to  establish  the  prior  uncertainty  regions  for  a,b  and  E.  The 
values  chosen  here  are  for  illustrative  purposes  only. 

In  some  cases,  it  may  be  useful  to  consider  additional  observations  or  judgements  on 
yield  such  as  would  be  provided,  for  example,  by  CORRTEX  monitoring,  assumed  to  be 
errors-in-variable  surrogates  for  yields  through  Equation  (2.1),  i.e. 

W,  =  w,  +  (,,  (2.7) 
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when*  the  errors  e,  are  assumed  independent  with  differing  but  known  variances  af.  Note 
that  assuming  the  variances  were  equal  over  calibrations  would  allow  estimating  them  as 
parameters  via  maximum  likelihood.  Consultations  with  scientific  sources  indicated  that 
this  assumption  would  not  be  a  reasonable  one  for  measured  hydrodynamic  yields.  It  seems 
natural  to  model  the  CORRTEX  observation  as  above  but  the  errors  in  variable  structure 
introduces  a  problem  in  the  updating  of  the  slope  and  intercept  matrix.  Because  of  the 
empirical  Bayes  nature  of  the  proposed  estimators  for  slope  and  intercept,  there  is  no  term 
in  the  equation  for  their  variances  that  involves  the  error  variance  cr2  of  the  errors  in  variable 
measurement  (see  Shumway,  1991).  Therefore,  we  have  taken,  as  an  alternate  approach,  the 
point  of  view  that  the  CORRTEX  yield  is  simply  another  observation  on  vector  magnitude 
with  mean  intercept  a0  =  0  and  mean  slope  b0  —  1  with  small  prior  variances.  Then,  the 
CORRTEX  measurement  fits  naturally  into  (2.1)  as  another  component,  of  magnitude'. 

The  incorporation  of  the  CORRTEX  yield  into  the  basic  model  means  that  we  must 
introduce  more  flexibility  into  the  priors  than  is  available  in  those  defined  in  Section  2.  For 
example,  the  priors  there  imply  that  the  joint  uncertainty  in  the  slope  and  intercept  vector 
is  proportional  to  the  uncertainty  in  the  error  in  the  magnitude  observation.  It  seems 
desirable  to  extend  the  possibilities  so  that  we  may  have,  for  example,  large  uncertainties 
on  the  CORRTEX  error  at  the  same  time  we  have  small  uncertainties  about  the  values 
a  =  0  and  6=1  for  CORRTEX.  What  we  give  up  with  this  formulation  is  the  ability 
to  be  uncertain  about  the  error  covariance  matrix  E  which  can’t  be  easily  integrated  out. 
Hence,  we  will  assume  that  the  covariance  matrix  of  each  error  is  known,  i.e. 

coe(ej)  =  E,  (2.8) 

In  the  approach  given  here,  we  assume  that  the  intercept  and  slope  vectors  a  and  b 
have  a  different  more  general  joint  normal  distribution  in  that  the  2p  X  1  stacked  vector 
3  =  (a',b')'  has  a  normal  distribution  with  mean  3o  =  (a^b^)'  and  known  2 p  x  2 p 
covariance  matrix 


where  the  covariances  of  the  components  a  and  b  appear  in  the  p  x  p  blocks.  The  form 
for  the  prior  density  function  is 


*2(0)  oc  exp{-l-(P  -  30)'T.-l(3  -  do)}- 


(2.9) 
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The  problems  of  interest  for  the  model  determined  by  (2.1)  are  again  calibration  and 
yield  estimation.  By  calibration  is  meant  the  estimation  of  the  slope  and  intercept  vectors 
b  and  a  and  the  covariance  matrix  £  in  (2.1)  from  a  small  data  set  consisting  of  paired 
vector  magnitudes  and  announced  yields  or  yield  surrogates.  For  yield  estimation,  we 
seek  the  best  estimator  for  log  yields  w  for  a  collection  of  observed  magnitudes.  The  next 
section  covers  calibration. 

3.  Calibration 

In  this  section,  we  discuss  first  in  3.1  the  calibration  of  the  regression  relation  (2.1) 
using  conventional  multivariate  regression.  The  classical  approach  has  the  advantage  that 
no  assumptions  are  made  about  a,  b  and  £  except  that  they  are  fixed  and  unknown.  The 
disadvantage  is  that  sample  sizes  are  almost  always  small,  leading  to  considerable  uncer¬ 
tainty  in  the  estimated  parameter  values.  In  section  3.2  we  consider  a  Bayesian  version 
of  regression  that  allows  incorporating  prior  assumptions  about  the  slope,  intercept  and 
yield-adjusted  magnitude  covariance  matrix  £.  If  input  is  allowed  as  to  the  probability 
distributions  for  the  parameters,  uncertainties  can  be  tightened  up  but  there  is  the  chance 
that  bad  initial  assumptions  will  contaminate  results.  The  possibility  of  interpreting  the 
yields  as  errors  in  variables  observations  is  considered  in  Section  4.3  under  the  prior  in¬ 
formation  assumed  in  the  Bayesian  approach.  In  this  case,  we  derive  maximum  likelihood 
estimators  for  the  unknown  fixed  yields  in  the  calibration  sample. 
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3.1  Classical  Regression 


Under  the  assumption  that  the  errors  in  (2.1)  are  independent  zero  mean  normal 
variables  with  p  X  p  covariance  matrix  E,  (see  Anderson,  1984)  the  maximum  likelihood 
estimator  for  the  slope-intercept  matrix  is 

B  =  YX’C~1  (3.1) 

where 


C  =  XX' 

Noting  the  partitioned  form  in  (2.2)  and  (2.3)  enables  writing  (3.1)  as 


h  =  E,(u,»  ~  «’)(m«  ~  ™) 

£,(  »’«  -  M')2 

and 


(3.2) 


(3.3) 


a  =  m  —  bte  (3.4) 

where  w  and  rh  are  the  mean  log-yield  and  magnitudes  respectively.  In  particular,  (3.3) 
and  (3.4)  show  the  analogy  with  the  case  where  p  =  1.  They  also  emphasize  that  the 
results  depend  on  a  collection  of  single  log-yield  variables,  u>,,  on  the  righthand  sides.  This 
simplifies  the  usual  multivariate  arguments  used  to  invert  the  predictive  distribution  for  a 
future  unknown  yield. 

The  covariance  matrix  S  can  be  taken  as  known  if  sufficient  information  is  available 
as  to  the  values  of  the  magnitude  variances  and  correlations.  One  may  also  compute  the 
maximum  likelihood  estimator 

t  =  (N-  2 )~l(Y  -  BX)(Y  -  BXy.  (3.5) 

Note  that  the  above  estimators  will  be  well  defined  only  when  there  are  enough  observations 
A'  >  p  +  2  to  make  the  resulting  normal  distributions  non-singular. 


3.2  Bayesian  Regression 


Now,  given  that  the  above  prior  information  can  be  specified,  it  is  natural  that  the 
estimates  for  the  slopes,  intercepts  and  error  covariance  may  change  from  those  specified 
by  the  maximum  likelihood  estimators  (3.1)  and  (3.5)  respectively.  To  handle  this,  simply 
take  the  joint  density  of  V,  D  and  E  and  derive  the  joint  conditional  (on  }')  marginal 
densities  of  B  and  E.  The  conditional  density  of  the  data  Y'  is  the  usual  likelihood 

L(Y"|B,E)oc  |E|-5*vexp{-^r  E_1(F  -  BX)(Y  -  BX)'}.  (3.6) 

Taking  the  posterior  means  of  the  density 

P(5,S|r)«I(F|B,S)7r(B|S)7r(S) 
leads  to  Bayesian  estimators  proportional  to 

Bk  =  (BC  +  B0D0)(C  +  Dq)  1  (3.7) 

where  C  is  defined  in  (3.2),  Bn  =  (a,v,b;v)  and 

v  v  =  (m  +  jV  +  1  -  p)~:  (mE0  +  (N  -  2)E  +  (B  -  Bo)(C~l  +  D~x  )~1{B  -  B0)'^ .  (3.8) 

The  subscript  N  differentiates  the  Bayes  estimators  from  the  maximum  likelihood  estima¬ 
tors  B  and  serves  as  a  reminder  that  the  Bayes  estimators  are  based  on  a  sample  of  size 
.V.  These  estimators  are  seen  to  combine  the  prior  panel  information  and  the  data  to  form 
a  combined  estimator  for  the  intercepts,  slopes  and  error  covariance  matrices. 

Using  the  alternate  prior  discussed  at  the  end  of  Section  2  with  the  prior  covariance 
on  the  intercept  slope  vector  ft  =  (a'b' )'  given  as  it  is  convenient  to  write  the  model 
(2.1 )  in  the  form 


where 


m,  =  X,ft  +  e, 


(3.9) 


Xi  =  (Ip,wiIp)  (3.10) 

is  the  ]>  x  2 p  partitioned  matrix  and  Ip  denotes  a  p  x  p  identity  matrix.  Under  the  assump¬ 
tions  given  above,  it  is  simple  to  show  that  for  I'  =  (mi, ... ,  m,v),  the  Bayesian  estimator 
for  the  slope  and  intercept  vector,  say  /?/v  =  E(ft\Y)  is 
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(3.11; 


whore 


/  N 
'«=! 


m,  + 


/  N  v  -> 

S«=  EA'.'E."'^.+S;’)  (3.12) 

'«  =  1  ' 

is  the  conditional  covariance  matrix. 

It  is  clear  that  there  are  two  routes  that  one  may  take  towards  estimating  the  slope 
from  an  uncalibrated  test  site.  The  first  is  to  insist  on  fully  weighting  the  available  obser¬ 
vations  and  use  the  ordinary  multivariate  regression  estimators  given  in  (3.1)  and  (3.5)  or 
(3.11)  when  the  covariance  matrices  of  the  errors  are  known.  If  a  reasonable  number  of 
calibration  events  are  furnished  and  if  one  believes  the  yields  for  these  events,  it  is  possible 
that  the  slope  and  intercept  estimators  will  be  reasonably  good;  the  variances  of  the  esti¬ 
mators  can  be  read  from  the  estimated  covariance  matrix  £©C_1 .  However,  there  may  be 
compelling  reasons  for  incorporating  prior  information  into  the  estimation  procedure.  The 
first  is  that  there  may  be  considerable  geophysical  information  and  expert  opinion  that  can 
be  used  to  narrow  down  the  possible  values  for  the  slope,  intercept  and  error  variance.  The 
second  reason  is  that  the  observed  data  may  have  been  deliberately  modified  to  contribute 
to  a  distorted  magnitude-yield  picture.  The  Bayesian  solution  has  the  advantage  of  being 
able  to  weight  optimally  the  two  components  of  information  contributed  by  the  expert 
opinion  (prior  distribution)  and  the  data. 
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3.3  Empirical  Bayes  With  CORRTEX 


We  now  consider  deriving  maximum  likelihood  estimators  for  the  yields  uq, . . . .  (ro¬ 
under  the  model  given  at  the  ends  of  section  2  and  section  3.3  where  the  covariance  matrix 
of  the  errors  is  assumed  known.  Consider  the  joint  likelihood  of  the  N  magnitudes  and 
hydrodynamic  CORRTEX  yields,  where  we  assume  that  the  last  element  of  the  vector 
m,  contains  the  CORRTEX  yield  W,,  defined  as  an  observation  satisfying  (2.7).  The 
likelihood  is  written  as 


N 


N 


=  II  / P(mi\P’w*)*2{0)  dl3  (3.13) 

«=i  i=i’' 

where  p(rrq|/?,it>i)  denotes  the  density  (conditional  on  /?)  of  the  ith  magnitude  and  7r2(/3) 
denotes  the  prior  given  by  (2.9).  We  assume  for  convenience  that  all  X,  are  known.  It 
is  clear  from  (2.9)  and  (3.9)  that  are  jointly  normal  with  mean  X,0u  and 

covariance  matrix  determined  by 


cov(  m, ,  m j ) 


(XiZfiXf  +  Ei,  if  i  =  j 
{  XiWj,  if  ifj. 


The  joint  likelihood  (3.13)  is  seen  from  the  form  of  X,  in  (3.10)  to  be  a  highly  nonlinear 
and  intractable  function  of  the  unknown  log  yields  uq , . . . ,  tejv  and  we  consider  maximizing 
it  by  indirect  means. 

The  EM  algorithm,  see  Dempster,  Laird  and  Rubin  (1977)  offers  a  considerable  sim¬ 
plification  to  the  usual  nonlinear  optimization  approaches.  If  we  were  to  observe  the  vector 
0,  the  complete-data  log-likelihood  for  this  problem  would  be  of  the  form 


log  I'  oc  -log  |E,|  -  \{0  -  0o)'Z?(0  ~  flo) 

IV  N  N 

-  y  Y,  bg  lS'l  ~  E(m'  "  X>PyZ7'( m.  -  X, 3)  (3.14) 

i= 1  i=l 

In  order  to  maximize  (3.13),  the  EM  algorithm  defines  an  iterative  sequence  consisting  of 
maximizing  successively  E(logL'\Y).  For  this  to  work,  we  need  expressions  for  £'(.:1|}') 
and  cov{0\Y]  which  are  available  through  Equations  (3.11)  and  (3.12).  Then,  define  J\  = 
(a'yy,  b'^)'  and  the  partitioned  conditional  covariance  matrix 
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62Ni  =  blvS-^Ar.  (3.16) 

The  above  equations  suggest  the  following  iterative  procedure'  for  computing  the  max¬ 
imum  likelihood  estimators. 

1.  Set  the  initial  yield  estimates  at  the  CORRTEX  estimates  Wi,i  =  1, . . . ,  N 

2.  Evaluate  (3^  and  E/v  using  Equations  (3.11),  (3.12)  and  the  current  estimated  log- 
yields. 

3.  Adjust  the  log-yields  using  Equations  ( 3.15)-( 3.16). 

3.  Return  to  Step  2.  until  convergence. 

It  is  clear  that  the  approach  of  this  section  generates  empirical  Bayes  estimators  for 
the  slope-intercept  vector  ft  since  the  conditional  expectations  are  evaluated  at  estimated 
yields.  These  empirical  Bayes  estimators  are  functions  of  the  maximum  likelihood  yield 
estimators  and  will  have  variances  that  increase  as  the  variances  of  the  maximum  likelihood 
estimators.  The  problem  with  the  maximum  likelihood  estimators  is  that  there  are  as  many 
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unknown  yields  as  there  are  vector  observations  so  that  an  asymptotic  result  will  not  be 
available. 

This  completes  the  specification  of  the  method  for  estimating  the  slope  and  inter¬ 
cept  vectors  along  with  the  error  covariance  matrix  using  classical  and  Bayes  regression 
procedures.  In  the  next  section,  we  consider  the  problems  involved  in  estimating  a  future 
unknown  yield  using  the  new  magnitude  vector  and  the  calibrated  regressions  derived  here. 

4.  Yield  Estimation  and  Uncertainties 

We  concentrate  now  on  the  problem  of  inferring  the  yield  when  the  new  magnitude 
vector  m  is  given  as 


m  =  a  +  bu>  +  e 

That  is,  we  are  interested  in  solving  the  above  for  the  yield  10"’  given  an  observation  on  a 
new  magnitude.  Of  course,  it  is  quite  obvious  that  given  a  fixed  known  a  and  b  and  the 
covariance  matrix  £,  the  maximum  likelihood  estimator  for  the  log-yield 


w  = 


b'E-1(m  —  a) 

b'S-Jb 


is  called  for  because  it  is  unbiased  and  has  variance 


(4.1) 


2  _  1 
*  b'Z-'b 


(4.2) 


For  example,  in  the  most  trivial  case,  with  independent  magnitudes,  one  obtains  the 
weighted  estimator 


w  =  cj (mi  -  «i)  4-  c2(w2  -  a2) 


where 


t r; 


Ci  = 


s+s' 


often  referred  to  as  the  unified  yield  estimator. 
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Furthermore,  in  the  general  case  given  above,  there  is  a  95%  confidence  interval  of 
the  form  tv  ±  1.96<7„,.  Since  w  denotes  log-yield,  say  logy,  the  limits  expressed  in  terms  of 
yield  y  are  ( U  ~ 1 ,  V )  where 


U  =  lO1’96^  (4.3) 

is  the  so-called  uncertainty  factor.  It  is  convenient  because  the  lower  and  upper  limits  for 
yield  are  obtained  by  dividing  or  multiplying  the  estimated  yield  by  the  uncertainty  factor. 

The  difficulties  associated  with  using  the  simple  inverted  unified  estimator  above  when 
a  and  b  are  estimated  from  calibration  data  have  been  discussed  at  length  in  the  statistical 
literature  (for  example,  see  Brown,  1982,  Hoadley,  1970,  Oman,  1988,  Picard  and  Bryson. 
1989)  and  we  will  only  say  that  even  under  joint  normality  of  a  and  b,  the  problem  is 
intractable  for  the  direct  estimator  given  in  (4.1).  An  approach  that  focuses  on  the  use 
of  the  predictive  distribution  of  m  to  derive  confidence  intervals  in  the  classical  case  and 
prediction  intervals  in  the  Bayesian  case  has  been  given  by  Brown  (1982).  We  use  Brown’s 
method  when  a  calibration  set  is  available  and  no  prior  information  can  be  assumed  for  the 
slope-intercept  and  covariance  parameters.  This  corresponds  to  the  calibration  information 
available  from  the  regression  approach  described  in  Section  3.1.  This  method,  which  is 
essentially  a  generalization  of  that  given  by  Fieller  (1954)  is  reviewed  in  Section  4.1. 

Prior  information  from  expert  panels  or  other  sources  can  be  used  to  incorporate  dis¬ 
tributional  assumptions  about  the  uncertainty  of  the  slope-intercept  matrix  B  and  the 
covariance  matrix  S  as  in  Section  3.2.  Brown  (1982)  describes  a  fully  Bayesian  approach 
that  sets  a  prior  distribution  on  the  log-yield  iv  as  well  and  then  computes  probability 
intervals  from  the  posterior  distribution  of  w.  For  example.  Picard  and  Bryson  (1992) 
follow  this  approach  with  a  uniform  distribution  assumed  a-priori  for  log-yield.  The  pos¬ 
terior  distribution  of  log-yield  must  then  be  computed  by  numerical  integration  methods. 
An  alternate  approach,  that  retains  the  Bayesian  flavor,  is  to  retain  the  assumption  that 
yields  are  fixed  unknown  constants  and  to  carry  through  the  process  of  inverting  the  pre¬ 
dictive  distribution  of  magnitude  as  in  the  classical  Fieller  methodology.  The  intervals 
obtained  then  have  the  same  interpretation  as  do  the  classical  ones.  This  approach,  which 
is  essentially  new,  is  given  in  Section  4.2. 

Finally,  in  the  approach  incorporating  CORRTEX  measurements  given  in  Section  3.3. 
we  obtain  maximum  likelihood  estimators  for  the  slope-intercept  matrix  B  which  can  be 
used  as  inputs  to  estimate  yield  using  a  modified  version  of  Equation  (4.1).  We  discuss 
this  approach  in  Section  4.3. 
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4.1  Classical  Prediction  Intervals 


We  develop  a  classical  interval  for  a  new  observation  by  noting  that  under  the  as¬ 
sumptions  of  Section  3.1,  with  in  a  new  magnitude  vector  and  the  new  unknown  yield  w. 
the  distribution  of  the  residual 


m  —  a  -  bie 

is  multivariate  normal  with  mean  0  and  covariance  matrix  (1  +  q{ii>))^l  where 

q(w)  =  (c11  +  2 c12u>  4-  c22w2)  (4-4) 

with  c,J  denoting  the  ijth  element  of  C~l .  It  is  also  the  case  that  the  residual  is  independent 
of  ( N  —  2)S,  which  has  the  Wishart  distribution  with  N-2  degrees  of  freedom.  Then,  the 
quadratic  form  in  the  residuals  with  the  inverse  of  the  covariance  matrix  (1  +  </(ir))L  in 
the  middle  has  Hotellings  T 2  distribution  as  in  Anderson  (1984,  p.  163)  or  Brown  (1982). 
Hence,  the  resulting  inequality 


(m  -  a  -  bu>)'£  *(m  -  a  -  bw) 


<  A'Pi,v(o) 


(4-5) 


(1  +q(w)) 

involving  the  quadratic  form  in  the  residuals  can  be  inverted  to  obtain  a  100(1  —  o  f'Z. 
confidence  interval.  The  constant  is 


PiN- 2)  r 

hp,N{a)  —  /  »r  .  ,  fp.N—  p-  1  (Q  b 


(4.6) 


(iV-p-1) 

with  Fp,/v-p-i(a)  denoting  the  (1  — a)  critical  point  on  the  F  distribution  with  p  numerator 
and  N-p-1  denominator  degrees  of  freedom.  We  obtain  the  ( 1  —  a  )  lev«'l  confidence  interval 
as 

d 

(4.7) 


d  r 
-±L 

c 


where 


ce 


(4.8) 


ilefines  the  upper  ami  lower  limit  with 


c  =  b'E-1b-c22A'p,N(a), 


(4.9) 
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and 


(4.9) 

(4.10) 


r  =  b,E-,b-r22A'J,.w(«), 

d  —  —  a)  -t-  c12  Kpx(a) 

e  =  (m-a)'£_1(m-a)-(l  +  c11  )A'Pijv(a)  (4.11) 

The  confidence  intervals,  derived  by  Brown  (1982),  are  the  exact  analogues  of  the 
univariate  intervals  first  given  by  Fieller  (1954).  As  in  the  univariate  case,  the  intervals 
may  be  open,  but  this  only  occurs  when  the  slope  vectors  are  not  significantly  different  from 
zero  at  level  n.  Oman  (19S8)  has  proposed  an  alternative  procedure  based  on  geometric 
arguments  when  there  is  a  strong  chance  that  the  intervals  above  might  be  degenerate. 

It  should  be  noted  that  the  interval  in  (4.7)  is  for  log-yield:  the  interval  on  yield  would 
be  obtained  by  inverting  the  limits  given  by  (4.7).  By  analogy  with  (4.3),  we  define  the 
uncertainty  number 


F=10L.  (4.12) 

Strictly  speaking,  this  requires  using  d/c  as  the  estimator  for  log-yield  but  we  will  custom¬ 
arily  use  the  slight  modification  of  this  implied  by  Equation  (4.1)  which  does  not  depend 
on  the  a  level  chosen  for  the  confidence. 

4.2  Predictive  Bayes  Methods 

We  may  also  look  at  deriving  posterior  predictive  intervals  under  the  assumption  that 
only  prior  information  is  available  in  the  form  of  the  joint  normal  inverted  Wishart  priors 
on  D  and  £  in  Section  2.  Now,  denoting  the  prior  inverted  Wishart  distribution  on  £  as 
jt(£)  and  the  conditional  density  of  the  intercept  and  slope  matrix  B  by  7t(Z?|E)  as  before, 
we  may  write  the  posterior  predictive  density  of  m  as 

P(m)  =  y*p(m|P,£)7r1(B|£)7r(£)  dB  d£  (4.13) 

where  we  integrate  first  over  B  and  then  over  £.  The  result  leads  to  a  multivariate  t 

distribution  for  m  (see,  for  example,  Anderson,  1984,  p.  283).  This  implies  that  the 
quadratic  form  in  the  multivariate  t  density  (involving  the  residuals  (m  —  ao  —  b0u’)  and 
the  inverse  of  a  particular  covariance  matrix)  has  the  F  distribution.  The  resulting  interval 
is  exactly  in  the  form  of  (4.8)  with 

c  =  bgEf'bo  -  do2A'p,m(a),  (4.14) 
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d  =  boS,  ‘(m  -  a0)  +  dlQ2Kp,m(a) 

and 

e  =  (m  -  ao/E'^m  -  a0)  -  (1  4-  d”  )A'p>m(a), 

where  dff  is  the  ijth  element  of  Dq1, 

s  _  ^Sp 

1  (m  ~  p  +  1) 

and 

A  piTn  ( O' )  =  pFp  m  —  p+ 1  (tt  )  (4.18) 

The  uncertainty  number  is  computed  as  before  using  (4.12)  with  the  modified  definitions 
of  the  constants  c,  d  and  e. 

Finally,  we  look  at  the  case  where  both  data  and  prior  information  are  available.  This 
means  that  we  introduce  the  joint  likelihood  of  the  calibration  data  Y  in  the  equation  for 
the  posterior  predictive  density.  This  leads  to 

P(m\Y)  =  J  p(m|B,E)Z(K|B,S)7r,(B|E)7r(S)  dB  dE  (4.19) 

where  L(Y\B ,  E)  denotes  the  likelihood  as  the  usual  conditional  density  of  the  data  given 
the  parameters.  Performing  the  above  integration  leads  to  the  posterior  t-distribution  as 
before  and  the  interval  (4.8)  obtains  again  with  the  values 

c  =  bNYN  bN-42/Vm,N(o),  (4.20) 

d  =  (m  —  a  at)  -+•  djy  Ap,m,/v(Q)  (4.21) 

and 

e  =  (m  —  ajv)'E^‘(m  -  ajV)  -  (1  +  d}J  )Kp.„,.i\( a).  (4.22) 

where  Bn  —  (a/v,b n)  so  that  and  b/v  can  be  picked  out  of  the  Bayesian  regression 
matrix  computed  in  Equation  (3.9)  in  Section  3.2  and  Ejv  is  defined  in  (3.10).  In  these 
expressions  d'^  denotes  the  ijth  element  of  D^1  where 


(4.15) 

(4.16) 

(4.17) 


D  n  —  C  +  Dq 


(4.23) 


and  the  constant  is 


Ap,m,/v(o? )  —  pF'p.m  +  N-\■^  —  j>(®  )• 


(4.24) 


16 


4.3  Predictive  Yield  Estimation  With  CORRTEX 

In  Section  3.3,  we  derived  maximum  likelihood  estimators  for  the  yields  in  a  calibration 
sample  using  CORRTEX  yields  as  the  last  magnitude  observation.  This  gave  empirical 
Bayes  estimators  for  the  slope-intercept  vector  /?/v  and  its  covariance  matrix  E,y  as  in 
Equations  (3.11)  and  (3.12)  where  we  evaluate  those  equations  at  the  estimated  log-yields 
in  the  calibration  sample.  Suppose  we  denote  these  empirical  Bayes  estimators  by  and 
E.y  respectively.  We  want  to  be  able  to  develop  an  estimator  of  a  single  new  yield,  given 
an  observation  vector  of  magnitudes  and  a  CORRTEX  yield  of  the  form 


m=a  +  bin  4-  e  (4.25) 

where  e  is  assumed  to  have  covariance  matrix  E  that  is  known. 

The  posterior  predictive  distribution  under  the  current  model  that  uses  the  prior 
density  tt2 ( 3 )  is  difficult  to  invert  in  order  to  derive  a  confidence  interval  for  the  yield 
measurement  as  was  done  in  Section  4.2  under  the  prior  density  tti(I?|E).  However,  it 
would  be  natural  to  use  an  estimator  like  (4.1)  in  this  case  and  we  consider  using  (4.1) 
evaluated  at  the  Bayes  estimators  a  /v ,  b  y ,  leading  to 


w  = 


b'^E  1  ( m  —  ) 

b'^E-’b* 


(4.26) 


We  can  show  that,  conditional  on  the  calibration  data  Y,  E(ib\Y)  =  w  and  that 


where 


£[(u>  -  u>)2|F]  = 


q(w) 
6%  ’ 


6%  =  b'^E  'b/v  (4.27) 

and  q(w)  is  a  quadratic  form  in  w.  The  posterior  distribution  is  normal  and  we  can  solve 
the  quadratic  equation 


{ib  -  w)2t>%  _2 

q(w)  ~ 

for  u\  where  ra.  denotes  the  upper  100(1  —  |)  percentile  on  the  standard  normal  dis¬ 
tribution.  We  will  use  o  =  .025  in  agreement  with  the  standard  used  in  previous  yield 
calculations. 
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The  above  algebra  leads  to  a  posterior  confidence  interval  of  the  form  (4.7)  where 


c  =  l-A^b'NE-1EN6S-1b/v  (4.28) 

d  =  w  —  A^b^E^EjvafcE^b/*/,  (4.29) 


and 


e  =  w 2  -  K2Q{Sl  +  b'yvE-1E^flS-1bjv),  (4.30) 

with 

A0  =  (4.31) 

°N 

The  uncertainty  factors  are  computed  as  before,  using  Equation  (4.12). 

The  practical  application  of  this  approach  would  use  an  initial  set  of  CORRTEX  yields 
to  derive  yield  estimators  using  the  maximum  likelihood  procedure  of  Section  3.3.  This 
also  yields  the  empirical  Bayes  estimators  a«  and  bv  for  the  intercept  and  slope  vec¬ 
tors  along  with  the  estimated  covariance  matrix  E/y-  Then,  new  yields  can  be  estimated 
using  Equation  (4.26)  with  only  seismic  values  or  with  combined  new  seismic  and  COR¬ 
RTEX  measurements.  If  only  seismic  values  are  available,  simply  use  the  submatrices  and 
subvectors  without  the  CORRTEX  parts  for  the  estimation  procedure. 
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5.  Calibration  and  Estimation:  An  Example 


In  order  to  give  a  simple  example  illustrating  the  computations  involved  in  the  various 
techniques  developed  here,  we  consider  a  set  of  16  Semipalatinsk  events  for  which  yields 
were  published  in  Bochaiov  et  al  (1989)  and  Vergino  (1989).  For  this  set  of  16  events, 
we  obtained  the  network  averaged  body  wave  magnitudes  mj,  corrected  for  receiver  terms 
as  well  as  the  near-source  focusing  and  defocusing  effects  from  Jih  and  Shumway  (1991). 
The  second  component  that  may  be  useful  is  root  mean  squared  error  computed  from  Lg 
magnitudes  at  a  network  of  Soviet  stations  at  regional  distances  ( 1000-4000  km)  as  given 
by  Israelsson  (1991).  The  16  events,  along  with  dates  are  shown  in  Table  1  below. 

Table  1:  Magnitudes  and  Announced  Yields  From  16  Explosions 


Event  No. 

Date 

mb 

l9 

Announced  Yield 

1 

03/20/66 

5.853 

6.048 

100 

2 

05/07/66 

4.456 

4.889 

4 

3 

09/29/68 

5.641 

5.801 

60 

4 

07/23/69 

5.186 

5.427 

16 

5 

11/30/69 

5.945 

6.024 

125 

6 

12/28/69 

5.660 

5.707 

46 

7 

04/25/71 

5.826 

5.946 

90 

8 

06/06/71 

5.319 

5.369 

16 

9 

10/09/71 

5.136 

5.192 

12 

10 

10/21/71 

5.341 

5.479 

23 

11 

02/10/72 

5.289 

5.369 

16 

12 

03/28/72 

4.984 

4.982 

6 

13 

08/16/72 

4.921 

4.948 

8 

14 

09/02/72 

4.602 

4.659 

2 

15 

11/02/72 

6.160 

6.158 

165 

16 

12/06/72 

5.983 

6.111 

140 

The  magnitudes  are  plotted  against  the  log-yields,  w  in  Figure  1  and  we  note  the 
excellent  agreement  with  the  linear  model  (2.1).  The  exception  is  the  mb  value  for  the 
second  event  which  seems  low  compared  to  the  Lg  and  the  other  points.  The  sample 
computations  here  will  assume  that  the  magnitude  vector  and  yield  values  are  available 
for  the  first  six  events  (above  the  line  in  the  table)  and  these  six  events  form  the  calibration 
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set.  The  symbols  for  the  six  events  are  shaded  in  Figure  1  and  we  see  that  they  cover  a 
reasonable  range  of  yields. 

We  assume  also  that  prior  information  in  the  form  of  expert  opinion  puts  the  intercepts 
of  the  two  lines  at  4.4  and  the  slopes  at  .9.  The  uncertainties  are  defined  as  in  Section 
2  by  assigning  standard  deviations  of  .2  to  the  intercept  vectors.  .1  to  the  slope  vectors 
and  a  correlation  of  —.5  between  the  intercept  and  slope.  We  always  assume  that  the 
and  Lg  yield  adjusted  magnitudes  have  standard  deviations  .05  and  .03  respectively  with 
correlation  .3.  This  matrix  becomes  So  in  the  Bayesian  approach  defined  by  the  first  set  of 
priors  given  in  (2.5)  and  (2.6).  In  this  formulation,  the  announced  yields  play  the  part  of  the 
known  wt  needed  for  the  approach  leading  to  (3.7)  and  (3.8).  Figure  1  shows  the  regression 
line  at  the  expected  values  and  it  is  clear  that  it  lies  above  the  line  that  might  be  implied 
by  the  data.  Inferring  yields  from  this  line  will  lead  to  severe  underestimates.  Figure  2 
shows  the  uncertainty  region  for  the  intercept  and  slope;the  intercept  is  between  3.9  and 
4.9  approximately  95%  of  the  time  with  the  comparable  assumptions  on  slope  involving 
the  interval  .65  to  1.15.  The  assumed  correlation  (slope  increases  as  the  intercept  gets 
smaller)  of  -.5  restricts  this  variation  somewhat. 

In  the  case  where  we  wish  to  incorporate  CORRTEX  values,  it  seems  reasonable  to 
regard  the  observation  on  CORRTEX  yield  as  another  element  in  the  magnitude  vector 
with  do  =  0  and  b0  —  1.  Then,  the  prior  distribution  given  by  (2.9)  has  three  components 
and  we  may  assigned  the  prior  standard  deviations  of  intercept  and  slope  for  the  COR¬ 
RTEX  magnitude  to  be  a  small  number  (.03  in  the  present  example)  with  no  correlation 
between  the  intercept  and  slope.  The  matrix  E#  in  (2.9)  is  a  6  x  6  matrix. 

Table  2  below  shows  the  estimated  intercepts  and  slopes  for  various  methods.  The 
first  two  lines  are  the  result  of  simple  multivariate  regressions  not  involving  the  use  of 
prior  information.  This  is  the  methodology  of  Section  3.1  and  we  see  that  the  calibration 
regression  on  the  first  six  events  differs  somewhat  from  the  regression  on  all  16  events.  In 
particular,  the  low  mj  value  pulls  the  intercept  down  to  3.92  and  increases  the  slope  to 
.98.  The  regression  on  the  full  set  of  16  events  gives  a  more  representative  line.  This  shows 
how  far  wrong  one  can  go  when  there  is  a  presumed  outlier  and  we  depend  on  the  small 
calibration  sample. 
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Table  2:  Estimated  Slopes  and  Intercepts  for  Semipalatinsk 


Mfc 

Intercepts 

Lg 

Mb 

Slopes 

Lg 

Data  (16  events) 

4.22 

4.34 

.83 

.83 

Data  (6  events) 

3.92 

4.45 

.98 

.77 

Bayes  Prior 

4.40 

4.40 

.90 

.90 

Bayes- Data  (6  events) 

3.97 

4.43 

.96 

.78 

Bayes- CORRTEX  (6  events) 

3.97 

4.45 

.96 

.78 

The  pure  Bayes  approach  calibrates  with  the  assumed  prior  distribution  and  these 
prior  assumptions  are  reproduced  in  row  3  of  Table  2.  Combining  the  prior  assumptions 
with  the  first  6  calibration  events  again  pulls  the  mj  intercept  down  and  the  slope  up  but 
the  Lg  slope  and  intercept  are  quite  close  to  what  obtains  for  the  entire  16  points.  All 
three  methods  give  nearly  the  same  answer  when  applied  to  the  calibration  set. 

For  the  CORRTEX  based  estimator,  we  use  the  empirical  Bayes  estimators  from 
Equation  (3.11),  with  yields  evaluated  by  the  maximum  likelihood  procedure  of  Section 
3.3.  The  uncertainties  are  evaluated  using  Equation  (3.12)  and  lead  to  the  estimated  stan¬ 
dard  deviations  .06,  .04,  .009,  .04,  .02,  .009  for  (01,02,03,61,62,63)  with  correlations  —.95 
between  cij ,  6|  and  02,62  and  .27  between  Oi,o2  and  between  6^ 62.  The  pairs  oi,62  and 
o2,6i  were  negatively  correlated  (  —  .26)  and  all  other  posterior  correlations  were  zero. 

The  result  of  applying  the  methodology  of  Section  3  to  get  yield  estimators  for  Events 
7-16  are  shown  in  Table  3.  One  can  compare  the  estimated  yields  with  the  announced 
yields  to  estimate  how  well  each  method  works.  The  data-driven  estimator  is  discussed  in 
Section  4.1.  The  estimator  is  (4.1)  evaluated  at  a,  b  and  the  endpoints  of  the  confidence 
interval  determine  the  uncertainty  number  using  (4.8)-(4.12).  The  pure  Bayes  estimator 
involves  only  prior  information  and  the  interval  defined  by  (4.14)-  (4.18)  with  the  estimator 
computed  by  evaluating  (4.1)  at  ao.  bo.  Combining  data  arid  the  prior  information  gives 
the  Bayes-data  estimator  by  evaluating  (4.1)  at  b/v.a/v  and  using  the  confidence  intervals 
defined  in  (4.20)-(4.24). 
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Table  3:  Yield  Estimation  Results  for  Semipalatinsk 
(6  calibration  events) 


Event 

Yield 

Data  Driven 

Pure  Bayes 

Bayes-Data 

CORRTEX 

7 

90*(2.01) 

8 

16 

9 

12 

11(1.86) 

7(2.28) 

10(1.63) 

110.28) 

10 

23 

23(1.94) 

15(2.24) 

23(1.59) 

23(1.31) 

11 

16 

18(1.89) 

12(2.27) 

17(1.60) 

17(1.30) 

12 

6 

6(1.75) 

4(2.33) 

6(1.68) 

6(1.26) 

13 

8 

5(1.82) 

4(2.33) 

5(1.69) 

6(1.26) 

14 

2 

2(1.73) 

1(2.48) 

2(1.83) 

2(1.22) 

15 

165 

185(2.12) 

98(2.56) 

170(1.66) 

165(1.39) 

16 

140 

145(2.07) 

81(2.47) 

140(1.65) 

135(1.38) 

*  Uncertainty  factors  shown  in  (  ) 


The  results  make  it  clear  that  the  pure  Bayes  apor  ,a<\i  goes  wrong  because  of  the 
choice  of  a  slightly  erroneous  prior;  all  yield?  '«n  underestimated  although  the  uncertainty 
number  is  large  enough  to  put  all  of  the  announced  yields  in  the  predicted  uncertainty 
region.  The  data-driven  regression  estimator  Jo; .3  reasonably  well  except  for  overestimat¬ 
ing  the  two  highest  yields;  again  the  uncertainty  intervals  include  the  true  values.  Both 
the  Bayes-data  and  CORRTEX  methods  do  quite  well.  The  Bayes-CORRTEX  estimator 
should  do  well  since  there  is  an  extra  observation  in  the  vector.  For  the  particular  example 
given  here,  the  Bayes-CORRTEX  method  assigned  weighting  coefficients  .21,. 60,  .33  to 
mk,Lg  and  the  CORRTEX  yield  respectively. 
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6.  Discussion 


There  are  many  plausible  approaches  to  calibration  and  estimation  in  the  multivariate  case. 
One  must  not  only  be  prepared  to  work  with  calibration  data  when  it  is  available  but  must 
be  willing  to  develop  a  yield  estimator  when  no  calibration  data  is  provided.  It  is  the  second 
of  these  two  situations  that  may  pertain  more  often  in  the  current  testing  environment 
where  monitoring  worldwide  proliferation  is  of  primary  interest.  In  such  an  environment, 
the  use  of  prior  information  provided  by  expert  panels  and  by  seismic  studies  of  different 
regions  of  the  earth  becomes  paramount.  The  emphasis  shifts  from  an  approach  that  is 
heavily  weighted  towards  calibration  from  samples,  possibly  involving  CORRTEX  yields, 
to  one  that  is  oriented  towards  formulating  more  accurate  prior  distributions.  Because  of 
the  emphasis  up  to  this  point  on  the  monitoring  of  a  threshold  test  ban  treat  (TTBT),  the 
research  summarized  in  this  report  concentrates  on  the  sampled  data  approach. 

However,  it  should  be  pointed  out  that  the  predictive  Bayes  method  discussed  at  the 
beginning  of  Section  4.2  is  an  example  of  an  approach  that  uses  only  prior  information  and 
needs  no  calibration  sample.  The  specification  of  the  prior  information  involves  expressing 
the  joint  uncertainty  in  the  slope-intercept  vectors  in  terms  of  a  multivariate  normal  distri¬ 
bution.  The  mean  value  of  this  multivariate  normal  (ao,bo)  contains  the  expected  slopes 
and  intercepts  whereas  the  uncertainty  is  specified  through  a  covariance  matrix.  Geomet¬ 
rically  this  involves  being  able  to  intuitively  choose  elliptical  uncertainty  contours  of  the 
form  given  in  Figure  3  that  express  our  prior  information  about  the  slope  and  intercept  . 

The  slope-intercept  prior  is  proportional  to  the  yield-  adjusted  magnitude  covariance 
matrix  which  may  also  be  subject  to  some  uncertainty.  This  is  expressed  through  the 
inverted  Wishart  distribution  with  a  specified  prior  covariance  matrix  given  by  Ed-  The 
marginal  distributions  behave  like  inverted  chi-  squared  distributions  with  the  uncertainty 
specified  by  the  degree  of  freedom  parameter  in.  Although  the  densities  are  restricted  in 
some  sense,  they  tend  to  emulate  what  one  might  regard  as  reasonable  prior  specifications 
for  uncertainty  about  the  parameter  E. 

If  the  above  prior  densities  are  not  flexible  enough,  the  approach  involving  the  second 
set  of  priors  given  by  Equation  (2.9)  can  be  applied.  This  does  impose  the  restriction 
that  we  must  be  willing  to  fix  the  yield  adjusted  magnitude  covariance  matrix  at  some 
known  value.  This  set  of  priors  seems  to  be  the  only  realistic  way  that  one  can  incorporate 
the  errors  in  variable  observation  on  CORRTEX  into  the  magnitude  vector.  The  posterior 
variances  of  the  slope-intercept  vector  will  be  weighted  properly  by  the  CORRTEX  variance 
in  this  particular  configuration. 
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Unfortunately,  any  Bayes  approach  involves  the  risk  of  using  an  incorrect  prior  so  that 
the  estimated  posterior  yield  might  be  quite  far  off.  However  the  uncertainty  factors  will 
increase  in  proportion  to  the  variance  covariance  structures  of  the  slope-intercept  vectors 
and  the  yield-adjusted  magnitude  vector.  Therefore,  the  overall  statements  including  the 
uncertainty  factor  will  include  the  true  yield  95%  of  the  time  if  the  uncertainty  of  the  prior 
intercepts  and  slopes  in  assigned  in  a  realistic  fashion.  This  was  seen  in  the  example  shown 
in  Table  3. 
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Priors  for  SD  (m  =  10,40) 


Std.  Dev. 


Figure  2:  Prior  Probability  Distributions  for  Standard  Deviation  for  a0  =  .05  (m=10  and 
40  degrees  of  freedom). 
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